Accelerometer-based prediction of skeletal mechanical loading during walking in normal weight to severely obese subjects:
Data analysis report
Data analysis report
Lucas Veras
2019-07-18
1 Background
Considering the prevalence of obesity (1), its associated fracture risk [(2)] and that these patients are especially prone to develop loading associated musculoskeletal injuries [(5)], there is a need to develop ground reaction force (GRF) prediction equations that are accurate for overweight and obese subjects in order to precisely determine and monitor exercise-associated mechanical loading in these patients. Therefore, the purpose of this study was to develop prediction equations based on accelerometer data able to accurately predict GRF on a broad range of body masses, from normal weight to severely obese subjects, setting thereby, a base for objective prescription and monitoring of exercise mechanical loads.
2 Load data
library(knitr)
library(broman)
library(tidyverse)
library(here)
library(nlme)
library(piecewiseSEM)
library(ez)
source(here("R", "cross_validate_mixed_model.R"))
source(here("R", "accuracy_indices.R"))
source(here("R", "get_BA_plot.R"))Initially, 3 files needed to be loaded, each containing data from one of the accelerometer placements utilized in the study: ankle, back and hip.
ankle <- read_csv(here("data", "ankle_sec.csv")) %>%
select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>%
filter(speed <= 6)
back <- read_csv(here("data", "back_sec.csv")) %>%
select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>%
filter(speed <= 6)
hip <- read_csv(here("data", "hip_sec.csv")) %>%
select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>%
filter(speed <= 6)A piece of a data frame used in all subsequent analysis can be seen below.
## # A tibble: 299 x 11
## ID speed body_mass height BMI sex age pVGRF_N pVACC_g pRGRF_N
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 3 2 85.6 171 29.3 M 50 1003. 2.21 1038.
## 2 3 3 85.6 171 29.3 M 50 1114. 2.84 1132.
## 3 3 4 85.6 171 29.3 M 50 1220. 3.64 1243.
## 4 3 5 85.6 171 29.3 M 50 1212. 3.78 1250.
## 5 3 6 85.6 171 29.3 M 50 1237. 3.76 1285.
## 6 4 2 61.1 162 23.3 F 45 670. 1.81 673.
## 7 4 3 61.1 162 23.3 F 45 694. 2.20 710.
## 8 4 4 61.1 162 23.3 F 45 743. 2.74 756.
## 9 4 5 61.1 162 23.3 F 45 803. 3.17 819.
## 10 4 6 61.1 162 23.3 F 45 872. 3.66 878.
## # … with 289 more rows, and 1 more variable: pRACC_g <dbl>3 Ground reaction force X acceleration relationship
We first determined the relationship between pACC and pGRF obtained during the incremental walking speeds used in our experimental protocol. This was performed in order to verify the consistency of this relationship for increasing pACC values and the ability of the accelerometer to discriminate differences in pGRF between subjects in different BMI classes. A scatterplot with these relationships for all three accelerometer placements tested and for both the resultant and its vertical component is depicted in Figure 1. Generally, as expected, pACC recorded by accelerometers in all placements, showed a linear increase as a function of pGRF increases. Also, the recorded accelerations were shown to be able to provide a good discrimination between different BMI classes, as for the same registered pACC the pGRF tended to be consistently higher for subjects in higher BMI classes.
Code to generate these plots can be seen here.
4 Linear mixed models
Linear mixed models (LMM) were used to estimate peak resultant ground reaction force (pRGRF) and peak vertical ground reaction force (pVGRF). Distinct LMMs were developed with data from ankle, back and hip accelerometer placement (ankle, back and hip data frames) using the lme() function of the nlme package. Covariance structure used was an autoregressive process of order 1 (correlation = corAR1) and maximum likelihood method was used for estimating parameters (method = "ML"). Predictors tested on pRGRF models were body mass and peak resultant acceleration (pRACC), while on pVGRF were body mass and vertical acceleration (pVACC). All of them were tested as fixed effects and have shown to be significant predictors (e.g.: fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass). Both random intercept and slopes were tested, but only the random intercept inclusion has showed models improvement (random = ~ 1 | ID). Linear, quadratic and cubic polynomial simulations were also tested, whereas the last one did not contribute significantly to the models. Final models were chosen according to -2 log-likelihood statistics. Traditional coefficient of determination (R2) was represented by conditional R2 (6) computed with rsquared() function of the piecewiseSEM package, that estimates the variance explained by the whole model (7).
Code to build the LMMs, as well as their summary and R2, can be found below.
# For resultant peak ground reaction force
ankle_res_LMM <- lme(
fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = ankle
)
r2_ankle_res_LMM <- rsquared(ankle_res_LMM)
summary(ankle_res_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: ankle
## AIC BIC logLik
## 3188.582 3218.186 -1586.291
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02697863 86.14005
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8748119
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) 188.55057 62.95028 232 2.995230 0.0030
## pRACC_g -2.04012 18.97302 232 -0.107527 0.9145
## I(pRACC_g^2) -6.97317 2.21342 232 -3.150405 0.0018
## body_mass 6.56845 0.70675 62 9.293826 0.0000
## pRACC_g:body_mass 1.88108 0.15506 232 12.130984 0.0000
## Correlation:
## (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g -0.700
## I(pRACC_g^2) 0.247 -0.700
## body_mass -0.898 0.426 0.130
## pRACC_g:body_mass 0.668 -0.562 -0.170 -0.766
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.59736240 -0.63264008 0.03559703 0.69274354 2.95031219
##
## Number of Observations: 299
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pRGRF_N gaussian identity none 0.9297268 0.9297268back_res_LMM <- lme(
fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back
)
r2_back_res_LMM <- rsquared(back_res_LMM)
summary(back_res_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back
## AIC BIC logLik
## 3185.997 3215.628 -1584.999
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.0184402 79.53301
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8530994
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -698.7031 123.44832 233 -5.659884 0e+00
## pRACC_g 1047.5129 122.40109 233 8.558036 0e+00
## I(pRACC_g^2) -345.2605 36.10073 233 -9.563809 0e+00
## body_mass 3.8294 1.06086 62 3.609729 6e-04
## pRACC_g:body_mass 6.0219 0.61893 233 9.729449 0e+00
## Correlation:
## (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g -0.896
## I(pRACC_g^2) 0.648 -0.892
## body_mass -0.636 0.268 0.143
## pRACC_g:body_mass 0.565 -0.281 -0.167 -0.922
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.0955967 -0.6001009 -0.0930397 0.5792458 3.1083008
##
## Number of Observations: 300
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pRGRF_N gaussian identity none 0.9361655 0.9361655hip_res_LMM <- lme(
fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip
)
r2_hip_res_LMM <- rsquared(hip_res_LMM)
summary(hip_res_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip
## AIC BIC logLik
## 3091.034 3120.637 -1537.517
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02080893 71.4172
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8662087
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -300.9909 89.71869 232 -3.354829 9e-04
## pRACC_g 522.6850 71.61994 232 7.298037 0e+00
## I(pRACC_g^2) -171.5606 16.91100 232 -10.144910 0e+00
## body_mass 3.9596 0.82234 62 4.814986 0e+00
## pRACC_g:body_mass 5.3671 0.41930 232 12.799992 0e+00
## Correlation:
## (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g -0.878
## I(pRACC_g^2) 0.588 -0.854
## body_mass -0.766 0.422 0.030
## pRACC_g:body_mass 0.678 -0.471 -0.038 -0.891
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.54273847 -0.59285086 0.05857753 0.70507501 2.51256437
##
## Number of Observations: 299
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pRGRF_N gaussian identity none 0.9488446 0.9488446# For vertical peak ground reaction force
ankle_vert_LMM <- lme(
fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = ankle
)
r2_ankle_vert_LMM <- rsquared(ankle_vert_LMM)
summary(ankle_vert_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: ankle
## AIC BIC logLik
## 3190.571 3220.175 -1587.286
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02781982 91.52233
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8928218
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) 173.06487 77.90449 232 2.221500 0.0273
## pVACC_g 17.98949 39.07008 232 0.460442 0.6456
## I(pVACC_g^2) -24.86971 6.71678 232 -3.702623 0.0003
## body_mass 5.35461 0.82437 62 6.495401 0.0000
## pVACC_g:body_mass 3.19203 0.27580 232 11.573831 0.0000
## Correlation:
## (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g -0.767
## I(pVACC_g^2) 0.393 -0.785
## body_mass -0.856 0.414 0.081
## pVACC_g:body_mass 0.680 -0.526 -0.092 -0.804
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.84130261 -0.57735799 0.02962816 0.55766690 2.61025073
##
## Number of Observations: 299
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pVGRF_N gaussian identity none 0.9177371 0.9177371back_vert_LMM <- lme(
fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back
)
r2_back_vert_LMM <- rsquared(back_vert_LMM)
summary(back_vert_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back
## AIC BIC logLik
## 3234.998 3264.628 -1609.499
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02508531 102.5201
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.9079104
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -776.8934 145.57197 233 -5.336834 0
## pVACC_g 1042.9052 146.96000 233 7.096524 0
## I(pVACC_g^2) -336.2115 44.82805 233 -7.500024 0
## body_mass 6.2132 1.22304 62 5.080120 0
## pVACC_g:body_mass 5.0805 0.72861 233 6.972936 0
## Correlation:
## (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g -0.889
## I(pVACC_g^2) 0.652 -0.894
## body_mass -0.688 0.327 0.066
## pVACC_g:body_mass 0.587 -0.341 -0.101 -0.891
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2287736 -0.6177471 -0.1216603 0.4585545 3.7585643
##
## Number of Observations: 300
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pVGRF_N gaussian identity none 0.8935621 0.8935621hip_vert_LMM <- lme(
fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip
)
r2_hip_vert_LMM <- rsquared(hip_vert_LMM)
summary(hip_vert_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip
## AIC BIC logLik
## 3153.414 3183.017 -1568.707
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02202366 84.32584
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8869442
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -435.7365 122.69794 232 -3.551294 5e-04
## pVACC_g 586.6627 112.66171 232 5.207295 0e+00
## I(pVACC_g^2) -188.9689 28.55314 232 -6.618146 0e+00
## body_mass 5.8047 0.97485 62 5.954511 0e+00
## pVACC_g:body_mass 4.9544 0.54067 232 9.163472 0e+00
## Correlation:
## (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g -0.910
## I(pVACC_g^2) 0.715 -0.907
## body_mass -0.759 0.471 -0.125
## pVACC_g:body_mass 0.688 -0.524 0.135 -0.888
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.95040403 -0.53289689 -0.04360304 0.45305786 4.09646433
##
## Number of Observations: 299
## Number of Groups: 64## Response family link method Marginal Conditional
## 1 pVGRF_N gaussian identity none 0.9263472 0.9263473In all models, conditional R2 values ranged from 0.89 to 0.95.
5 Validation analysis
Model validation was assessed by the leave-one-out cross-validation (LOOCV) method. This approach is recommended when a different sample is not available for cross-validation and it provides an insight on the model potential to predict outcomes in a new independent sample (8). For LOOCV each participant’s data was separated into a testing dataset (one at a time) with the remaining data being in the training dataset. New LMM, with the same outcomes and predictors as determined for the entire sample, were developed using the training dataset and then used to predict pGRF for the only participant in the testing dataset. This process was repeated for all participants (64 times). Data from testing dataset were used in the remaining statistical analysis.
LOOCV of the LMM was done with the cross_validate_mixed_model() function and is shown below.
# For resultant peak ground reaction force
fix_eff <- pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass
rand_eff <- ~ 1 | ID
# Ankle
LOOCV_ankle_res_LMM <- do.call(rbind, (lapply(unique(ankle$ID), cross_validate_mixed_model, df = ankle)))
# Back
LOOCV_back_res_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_res_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))
# For vertical peak ground reaction force
fix_eff <- pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass
rand_eff <- ~ 1 | ID
# Ankle
LOOCV_ankle_vert_LMM <- do.call(rbind, (lapply(unique(ankle$ID), cross_validate_mixed_model, df = ankle)))
# Back
LOOCV_back_vert_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_vert_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))5.1 Bland-Altman plots
Bland-Altman plots were used to examine the agreement between pGRF measured with force plates and those predicted through the regression equations. The difference of the actual and predicted pGRF was plotted against their mean. Bias was expressed as the mean of these differences and the limits of agreement were obtained using ±1.96 standard deviation of the mean between actual and predicted pGRF (9).
Bland-Altman plots can be seen below for all accelerometer placements. Panel A shows plots of pRGRF and panel B shows plots of pVGRF.
Code to generate these plots can be seen here.
One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.
# For resultant peak ground reaction force
# Ankle
LOOCV_ankle_res_LMM$diff <- LOOCV_ankle_res_LMM$pRGRF_N - LOOCV_ankle_res_LMM$pRGRF_N_predicted
LOOCV_ankle_res_LMM$mean <- (LOOCV_ankle_res_LMM$pRGRF_N + LOOCV_ankle_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_ankle_res_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_ankle_res_LMM$diff
## t = 0.37558, df = 298, p-value = 0.7075
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -9.266825 13.638154
## sample estimates:
## mean of x
## 2.185665# Back
LOOCV_back_res_LMM$diff <- LOOCV_back_res_LMM$pRGRF_N - LOOCV_back_res_LMM$pRGRF_N_predicted
LOOCV_back_res_LMM$mean <- (LOOCV_back_res_LMM$pRGRF_N + LOOCV_back_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_back_res_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_back_res_LMM$diff
## t = 0.44993, df = 299, p-value = 0.6531
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.724802 12.303988
## sample estimates:
## mean of x
## 2.289593# Hip
LOOCV_hip_res_LMM$diff <- LOOCV_hip_res_LMM$pRGRF_N - LOOCV_hip_res_LMM$pRGRF_N_predicted
LOOCV_hip_res_LMM$mean <- (LOOCV_hip_res_LMM$pRGRF_N + LOOCV_hip_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_hip_res_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_hip_res_LMM$diff
## t = 0.36236, df = 298, p-value = 0.7173
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.774785 11.284133
## sample estimates:
## mean of x
## 1.754674# For vertical peak ground reaction force
# Ankle
LOOCV_ankle_vert_LMM$diff <- LOOCV_ankle_vert_LMM$pVGRF_N - LOOCV_ankle_vert_LMM$pVGRF_N_predicted
LOOCV_ankle_vert_LMM$mean <- (LOOCV_ankle_vert_LMM$pVGRF_N + LOOCV_ankle_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_ankle_vert_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_ankle_vert_LMM$diff
## t = -0.77815, df = 298, p-value = 0.4371
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -17.18212 7.44448
## sample estimates:
## mean of x
## -4.868821# Back
LOOCV_back_vert_LMM$diff <- LOOCV_back_vert_LMM$pVGRF_N - LOOCV_back_vert_LMM$pVGRF_N_predicted
LOOCV_back_vert_LMM$mean <- (LOOCV_back_vert_LMM$pVGRF_N + LOOCV_back_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_back_vert_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_back_vert_LMM$diff
## t = -0.10563, df = 299, p-value = 0.9159
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -12.19819 10.95537
## sample estimates:
## mean of x
## -0.6214066# Hip
LOOCV_hip_vert_LMM$diff <- LOOCV_hip_vert_LMM$pVGRF_N - LOOCV_hip_vert_LMM$pVGRF_N_predicted
LOOCV_hip_vert_LMM$mean <- (LOOCV_hip_vert_LMM$pVGRF_N + LOOCV_hip_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_hip_vert_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_hip_vert_LMM$diff
## t = 0.0087888, df = 298, p-value = 0.993
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -9.956701 10.046032
## sample estimates:
## mean of x
## 0.04466578Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.
For lower back and hip placements, results showed a significant proportional bias (p < 0.05), however, with a low magnitude (highest R2 = 0.032). These results suggest that despite there is a trend for underestimation at increasingly higher pGRF values, the magnitude of this effect is neglectable. No proportional bias (p > 0.05) was detected for ankle derived equations.
# For resultant peak ground reaction force
# Ankle
ankle_res_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_ankle_res_LMM)
summary(ankle_res_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_ankle_res_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -338.28 -54.43 3.59 57.89 285.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.03516 22.46494 -1.471 0.142
## mean 0.03011 0.01855 1.623 0.106
##
## Residual standard error: 100.4 on 297 degrees of freedom
## Multiple R-squared: 0.00879, Adjusted R-squared: 0.005453
## F-statistic: 2.634 on 1 and 297 DF, p-value: 0.1057# Back
back_res_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_res_LMM)
summary(back_res_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_res_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -267.995 -55.071 -2.505 52.450 294.243
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -47.77144 19.69289 -2.426 0.01587 *
## mean 0.04282 0.01628 2.630 0.00899 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.28 on 298 degrees of freedom
## Multiple R-squared: 0.02268, Adjusted R-squared: 0.0194
## F-statistic: 6.915 on 1 and 298 DF, p-value: 0.008991##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_res_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -268.721 -53.909 5.937 51.546 235.762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -40.38734 18.75117 -2.154 0.0321 *
## mean 0.03610 0.01552 2.325 0.0207 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 83.12 on 297 degrees of freedom
## Multiple R-squared: 0.01788, Adjusted R-squared: 0.01457
## F-statistic: 5.406 on 1 and 297 DF, p-value: 0.02074# For vertical peak ground reaction force
# Ankle
ankle_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_ankle_vert_LMM)
summary(ankle_vert_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_ankle_vert_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -418.05 -60.21 10.54 59.08 293.87
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -45.43603 24.52661 -1.853 0.0649 .
## mean 0.03523 0.02060 1.710 0.0883 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 107.8 on 297 degrees of freedom
## Multiple R-squared: 0.009752, Adjusted R-squared: 0.006418
## F-statistic: 2.925 on 1 and 297 DF, p-value: 0.08827# Back
back_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_vert_LMM)
summary(back_vert_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_vert_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -342.85 -64.14 -4.92 51.55 354.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -46.47386 23.16029 -2.007 0.0457 *
## mean 0.03998 0.01954 2.046 0.0416 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 101.4 on 298 degrees of freedom
## Multiple R-squared: 0.01386, Adjusted R-squared: 0.01055
## F-statistic: 4.187 on 1 and 298 DF, p-value: 0.04162# Hip
hip_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_vert_LMM)
summary(hip_vert_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_vert_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -255.386 -53.520 -0.846 49.427 295.833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -42.10744 19.93510 -2.112 0.0355 *
## mean 0.03683 0.01685 2.186 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.33 on 297 degrees of freedom
## Multiple R-squared: 0.01583, Adjusted R-squared: 0.01252
## F-statistic: 4.778 on 1 and 297 DF, p-value: 0.029615.2 Indices of accuracy
To evaluate models prediction accuracy, mean absolute error (MAE), mean absolute percent error (MAPE) and root mean square error (RMSE) were computed These indices were computed with the accuracy_indices() function.
# For resultant peak ground reaction force
# Ankle
ankle_res_accuracy <- accuracy_indices(LOOCV_ankle_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# Back
back_res_accuracy <- accuracy_indices(LOOCV_back_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# Hip
hip_res_accuracy <- accuracy_indices(LOOCV_hip_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# For vertical peak ground reaction force
# Ankle
ankle_vert_accuracy <- accuracy_indices(LOOCV_ankle_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")
# Back
back_vert_accuracy <- accuracy_indices(LOOCV_back_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")
# Hip
hip_vert_accuracy <- accuracy_indices(LOOCV_hip_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")Table below shows these indices for each accelerometer placement in both pRGRF and pVGRF.
accuracy <- data.frame(
acc_placement = c(rep("Ankle", 2), rep("Back", 2), rep("Hip", 2)),
condition = rep(c("pRGRF", "pVGRF"), 3),
MAE = c(
round(ankle_res_accuracy[[5]], 1), round(ankle_vert_accuracy[[5]], 1),
round(back_res_accuracy[[5]], 1), round(back_vert_accuracy[[5]], 1),
round(hip_res_accuracy[[5]], 1), round(hip_vert_accuracy[[5]], 1)
),
MAPE = c(
paste(round(ankle_res_accuracy[[12]], 1), ".0%", sep = ""), paste(round(ankle_vert_accuracy[[12]], 1), "%", sep = ""),
paste(round(back_res_accuracy[[12]], 1), "%", sep = ""), paste(round(back_vert_accuracy[[12]], 1), "%", sep = ""),
paste(round(hip_res_accuracy[[12]], 1), "%", sep = ""), paste(round(hip_vert_accuracy[[12]], 1), "%", sep = "")
),
RMSE = c(
round(ankle_res_accuracy[[7]], 1), round(ankle_vert_accuracy[[7]], 1),
round(back_res_accuracy[[7]], 1), round(back_vert_accuracy[[7]], 1),
round(hip_res_accuracy[[7]], 1), round(hip_vert_accuracy[[7]], 1)
)
)
kable(accuracy, col.names = c("Accelerometer Placement", "GRF", "MAE", "MAPE", "RMSE"), align = c("l", "l", "r", "r", "r"))| Accelerometer Placement | GRF | MAE | MAPE | RMSE |
|---|---|---|---|---|
| Ankle | pRGRF | 77.0 | 7.0% | 100.5 |
| Ankle | pVGRF | 81.0 | 7.5% | 108.1 |
| Back | pRGRF | 66.7 | 5.9% | 88.0 |
| Back | pVGRF | 74.9 | 6.5% | 101.7 |
| Hip | pRGRF | 65.5 | 5.9% | 83.6 |
| Hip | pVGRF | 65.6 | 5.9% | 87.7 |
All conditions showed a MAE from near 65 N to near 81 N, MAPE from near 6% to near 7% and RMSE from near 84 N to 84 N.
5.3 ANOVA
A series of repeated measures analysis of variance (ANOVA) were run to assess whether pGRF predicted from the regression equations were significantly different from those measured with FP. Walking speeds, accelerometer placements (ankle, back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis. These procedures were taken separately for resultant and its vertical component. If assumptions of sphericity were violated (p < 0.05), the conservative Greenhouse–Geisser correction factor would be applied to adjust the degrees of freedom. Post-hoc analyses would be conducted using pairwise comparisons with Holm’s test if a significant difference were observed among actual and predicted pGRF.
In order to run the repeated measures ANOVA, new data frames needed to be built, for both pRGRF and pVGRF, putting all the pGRF values in a single variable and grouping then by accelerometer placement or actual value in another variable. An example of code to tidy the data is shown below.
## Predicted pRGRF
ankle_res_pred <- LOOCV_ankle_res_LMM %>%
select(ID, speed, pRGRF_N_predicted) %>%
spread(key = speed, value = pRGRF_N_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = ankle
)
back_res_pred <- LOOCV_back_res_LMM %>%
select(ID, speed, pRGRF_N_predicted) %>%
spread(key = speed, value = pRGRF_N_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = back
)
hip_res_pred <- LOOCV_hip_res_LMM %>%
select(ID, speed, pRGRF_N_predicted) %>%
spread(key = speed, value = pRGRF_N_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = hip
)
res_pred_df <- ankle_res_pred %>%
full_join(back_res_pred, by = c("ID", "speed")) %>%
full_join(hip_res_pred, by = c("ID", "speed")) %>%
na.omit()
## Actual pRGRF
ankle_res_actual <- LOOCV_ankle_res_LMM %>%
select(ID, speed, pRGRF_N) %>%
spread(key = speed, value = pRGRF_N) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_ankle
)
back_res_actual <- LOOCV_back_res_LMM %>%
select(ID, speed, pRGRF_N) %>%
spread(key = speed, value = pRGRF_N) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_back
)
hip_res_actual <- LOOCV_hip_res_LMM %>%
select(ID, speed, pRGRF_N) %>%
spread(key = speed, value = pRGRF_N) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_hip
)
res_actual_df <- ankle_res_actual %>%
full_join(back_res_actual, by = c("ID", "speed")) %>%
full_join(hip_res_actual, by = c("ID", "speed")) %>%
na.omit()
res_actual_df <- res_actual_df %>%
mutate(actual = (actual_ankle + actual_back + actual_hip) / 3) %>%
select(ID, speed, actual)
## Merge predicted and actual data frames
res_ANOVA_df <- res_pred_df %>%
full_join(res_actual_df, by = c("ID", "speed")) %>%
gather(
ankle, back, hip, actual,
key = "group",
value = "pRGRF"
)
res_ANOVA_df$ID <- as.factor(res_ANOVA_df$ID)
res_ANOVA_df$speed <- as.factor(res_ANOVA_df$speed)
res_ANOVA_df$group <- as.factor(res_ANOVA_df$group)Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Summary statistics for pRGRF and pVGRF ANOVA can be seen below, respectively.
# For resultant peak ground reaction force
res_ANOVA <- ezANOVA(
data = res_ANOVA_df,
dv = pRGRF,
wid = ID,
within = .(speed, group),
detailed = TRUE,
type = 3
)
res_ANOVA## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 45 1.121949e+09 55569779.5 9.085464e+02 1.758136e-31
## 2 speed 4 180 1.012027e+07 439209.6 1.036890e+03 4.502391e-123
## 3 group 3 135 1.815789e+03 1484985.0 5.502448e-02 9.829444e-01
## 4 speed:group 12 540 2.892448e+04 508964.5 2.557352e+00 2.676332e-03
## p<.05 ges
## 1 * 9.508430e-01
## 2 * 1.485583e-01
## 3 3.130415e-05
## 4 * 4.984240e-04
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 speed 3.520032e-02 1.041884e-26 *
## 3 group 3.731595e-01 3.597677e-08 *
## 4 speed:group 3.580631e-06 1.091425e-64 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 speed 0.3944080 3.843045e-50 * 0.4063060 1.412612e-51
## 3 group 0.6892659 9.505875e-01 0.7229903 9.560913e-01
## 4 speed:group 0.2458901 5.887449e-02 0.2650287 5.417688e-02
## p[HF]<.05
## 2 *
## 3
## 4# For vertical peak ground reaction force
vert_ANOVA <- ezANOVA(
data = vert_ANOVA_df,
dv = pVGRF,
wid = ID,
within = .(speed, group),
detailed = TRUE,
type = 3
)
vert_ANOVA## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 45 1.082599e+09 52600265.0 926.1727721 1.163820e-31
## 2 speed 4 180 7.572969e+06 364333.1 935.3627935 3.185492e-119
## 3 group 3 135 1.658433e+04 2164731.4 0.3447516 7.929979e-01
## 4 speed:group 12 540 2.257064e+04 474638.4 2.1399001 1.338656e-02
## p<.05 ges
## 1 * 0.9511475606
## 2 * 0.1198692000
## 3 0.0002981690
## 4 * 0.0004057531
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 speed 4.088323e-02 2.251569e-25 *
## 3 group 6.550476e-01 2.393741e-03 *
## 4 speed:group 8.739178e-06 4.965719e-58 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 speed 0.4424565 3.134415e-54 * 0.4592394 3.449388e-56
## 3 group 0.7989465 7.473983e-01 0.8469999 7.594651e-01
## 4 speed:group 0.2757896 9.127102e-02 0.3002998 8.521643e-02
## p[HF]<.05
## 2 *
## 3
## 4For either resultant and its vertical component, repeated measures ANOVA demonstrated that actual and predicted pGRF increased significantly (p < 0.001) along with increments in walking speed.
To assess for any significant differences between actual and predicted pGRF in each speed, five separate repeated measures ANOVA were done, one for each walking speed.
# For resultant peak ground reaction force
# 2 km/h
res_ANOVA_s2 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 2),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 3 km/h
res_ANOVA_s3 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 3),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 4 km/h
res_ANOVA_s4 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 4),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 5 km/h
res_ANOVA_s5 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 5),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 6 km/h
res_ANOVA_s6 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 6),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# For vertical peak ground reaction force
# 2 km/h
vert_ANOVA_s2 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 2),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 3 km/h
vert_ANOVA_s3 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 3),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 4 km/h
vert_ANOVA_s4 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 4),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 5 km/h
vert_ANOVA_s5 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 5),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 6 km/h
vert_ANOVA_s6 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 6),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)No significant differences (p > 0.05) between actual and predicted pGRF in each speed were found.
These results can be seen in figure below.
6 Regressions comparison
Our equation was compared with a previously published reference equation (11), in which Neugebauer et al. used a similar approach for the prediction of pGRF. Comparison was performed using the equation to predict pVGRF from hip-worn accelerometers, as it was the only suitable with our data. The analysis was performed in three ways using: i) a subsample of participants with normal weight or overweight (BMI ≥ 18.5 and < 30 kg.m-2); ii) a subsample of obese participants (BMI ≥ 30 kg.m-2); iii) the whole sample.
New data frames needed to be built for each of the subsamples and the whole sample. Then, pVGRF was predicted using the reference equation. Code to build the data frames is shown below.
# Prepare data frames
non_obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>%
filter(BMI < 30)
obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>%
filter(BMI >= 30)
whole_sample_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted)
# Apply Neugebauer 2014 equation
non_obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(non_obese_df)) {
non_obese_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * non_obese_df$pVACC_g[i]) + (0.014 * non_obese_df$body_mass[i]))
}
obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(obese_df)) {
obese_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * obese_df$pVACC_g[i]) + (0.014 * obese_df$body_mass[i]))
}
whole_sample_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(whole_sample_df)) {
whole_sample_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * whole_sample_df$pVACC_g[i]) + (0.014 * whole_sample_df$body_mass[i]))
}Bland-Altman plots were used to confront the agreement between pVGRF measured with FP and those predicted using both regression equations.
Code to generate these plots can be seen here.
A lower dispersion around the bias value, as well as a lower percentage of values falling off the LoA can be appreciated in our equation. Additionally, while the bias in our equation tended to be zero, in the reference equation bias was always > 0 (p < 0.05), showing a consistent underestimation of pVGRF.
Also, to assess the prediction accuracy, MAE, MAPE and RMSE were calculated using pVGRF values predicted from both equations.
non_obese_our_accuracy <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_predicted")
non_obese_Neug_accuracy <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
obese_our_accuracy <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_predicted")
obese_Neug_accuracy <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
whole_sample_our_accuracy <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_predicted")
whole_sample_Neug_accuracy <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_Neugebauer")The accuracy indices from our equation were all substantially lower compared to the reference equation indices, with an overall MAPE approximately 3 times smaller for our equation. For both equations, the MAPE was lower in the obese subsample than in the non-obese subsample
7 R session info
## ─ Session info ──────────────────────────────────────────────────────────
## setting value
## version R version 3.5.2 (2018-12-20)
## os macOS Mojave 10.14.5
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Maceio
## date 2019-07-18
##
## ─ Packages ──────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 3.5.0)
## assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.0)
## backports 1.1.2 2017-12-13 [1] CRAN (R 3.5.0)
## broman * 0.68-2 2018-07-25 [1] CRAN (R 3.5.0)
## broom 0.4.5 2018-07-03 [1] CRAN (R 3.5.0)
## callr 2.0.4 2018-05-15 [1] CRAN (R 3.5.0)
## car 3.0-2 2018-08-23 [1] CRAN (R 3.5.0)
## carData 3.0-2 2018-09-30 [1] CRAN (R 3.5.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.5.0)
## cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.0)
## colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.0)
## cowplot * 0.9.3 2018-07-15 [1] CRAN (R 3.5.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.0)
## curl 3.2 2018-03-28 [1] CRAN (R 3.5.0)
## data.table 1.12.0 2019-01-13 [1] CRAN (R 3.5.2)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.0)
## devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.0)
## digest 0.6.19 2019-05-20 [1] CRAN (R 3.5.2)
## dplyr * 0.8.0.1 2019-02-15 [1] CRAN (R 3.5.2)
## evaluate 0.10.1 2017-06-24 [1] CRAN (R 3.5.0)
## ez * 4.4-0 2016-11-02 [1] CRAN (R 3.5.0)
## fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.0)
## forcats * 0.3.0 2018-02-19 [1] CRAN (R 3.5.0)
## foreign 0.8-71 2018-07-20 [1] CRAN (R 3.5.2)
## fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.0)
## ggplot2 * 3.1.0 2018-10-25 [1] CRAN (R 3.5.0)
## glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.0)
## gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.0)
## haven 2.0.0 2018-11-22 [1] CRAN (R 3.5.0)
## here * 0.1 2017-05-28 [1] CRAN (R 3.5.0)
## highr 0.7 2018-06-09 [1] CRAN (R 3.5.0)
## hms 0.4.2 2018-03-10 [1] CRAN (R 3.5.0)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.0)
## httr 1.3.1 2017-08-20 [1] CRAN (R 3.5.0)
## jsonlite 1.5 2017-06-01 [1] CRAN (R 3.5.0)
## knitr * 1.21 2018-12-10 [1] CRAN (R 3.5.0)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.5.2)
## lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.0)
## lme4 1.1-19 2018-11-10 [1] CRAN (R 3.5.0)
## lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.5.0)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.0)
## MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.2)
## Matrix 1.2-15 2018-11-01 [1] CRAN (R 3.5.2)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.0)
## mgcv 1.8-26 2018-11-21 [1] CRAN (R 3.5.2)
## minqa 1.2.4 2014-10-09 [1] CRAN (R 3.5.0)
## mnormt 1.5-5 2016-10-15 [1] CRAN (R 3.5.0)
## modelr 0.1.2 2018-05-11 [1] CRAN (R 3.5.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.0)
## nlme * 3.1-137 2018-04-07 [1] CRAN (R 3.5.2)
## nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.5.0)
## openxlsx 4.1.0.1 2019-05-28 [1] CRAN (R 3.5.2)
## pbkrtest 0.4-7 2017-03-15 [1] CRAN (R 3.5.0)
## piecewiseSEM * 2.0.2 2018-07-24 [1] CRAN (R 3.5.0)
## pillar 1.3.1 2018-12-15 [1] CRAN (R 3.5.0)
## pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.0)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.0)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.0)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.0)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.0)
## processx 3.1.0 2018-05-15 [1] CRAN (R 3.5.0)
## psych 1.8.4 2018-05-06 [1] CRAN (R 3.5.0)
## purrr * 0.3.0 2019-01-27 [1] CRAN (R 3.5.2)
## R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.2)
## Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.0)
## readr * 1.3.0 2018-12-11 [1] CRAN (R 3.5.0)
## readxl 1.1.0 2018-04-20 [1] CRAN (R 3.5.0)
## remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.0)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.0)
## rio 0.5.16 2018-11-26 [1] CRAN (R 3.5.0)
## rlang 0.3.1 2019-01-08 [1] CRAN (R 3.5.2)
## rmarkdown 1.10 2018-06-11 [1] CRAN (R 3.5.0)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.0)
## rstudioapi 0.7 2017-09-07 [1] CRAN (R 3.5.0)
## rvest 0.3.2 2016-06-17 [1] CRAN (R 3.5.0)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.0)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.0)
## stringi 1.3.1 2019-02-13 [1] CRAN (R 3.5.2)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.5.2)
## testthat 2.0.0 2017-12-13 [1] CRAN (R 3.5.0)
## tibble * 2.0.1 2019-01-12 [1] CRAN (R 3.5.2)
## tidyr * 0.8.2 2018-10-28 [1] CRAN (R 3.5.0)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.0)
## tidyverse * 1.2.1 2017-11-14 [1] CRAN (R 3.5.0)
## usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.0)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 3.5.0)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.0)
## xfun 0.4 2018-10-23 [1] CRAN (R 3.5.0)
## xml2 1.2.0 2018-01-24 [1] CRAN (R 3.5.0)
## yaml 2.1.19 2018-05-01 [1] CRAN (R 3.5.0)
## zip 1.0.0 2017-04-25 [1] CRAN (R 3.5.0)
##
## [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/libraryReferences
1. Guthold R, Stevens GA, Riley LM, Bull FC. Worldwide trends in insufficient physical activity from 2001 to 2016: a pooled analysis of 358 population-based surveys with 19 million participants. Lancet Glob Heal. 2018;6(10):e1077–86.
2. Li X, Gong X, Jiang W. Abdominal obesity and risk of hip fracture: a meta-analysis of prospective studies. Osteoporos Int. 2017;28(10):2747–57.
3. Hind K, Rudman H, Pearce M, Treadgold L, Birrell F. Obesity, bone density for weight and prevalent vertebral fracture at age 62 years: the Newcastle Thousand Families Study. J Clin Densitom. 2018;21(4):606.
4. Kaze AD, Rosen HN, Paik JM. A meta-analysis of the association between body mass index and risk of vertebral fracture. Osteoporos Int. 2018;29(1):31–9.
5. Viester L, Verhagen EA, Hengel KM, Koppes LL, Van Der Beek AJ, Bongers PM. The relation between body mass index and musculoskeletal symptoms in the working population. BMC Musculoskelet Disord. 2013;14(238).
6. Nakagawa S, Schielzeth H, O’Hara RB. A general and simple method for obtainingR2from generalized linear mixed-effects models. Methods in Ecology and Evolution. 2013;4(2):133–42.
7. Jaeger BC, Edwards LJ, Das K, Sen PK. An r2 statistic for fixed effects in the generalized linear mixed model. Journal of Applied Statistics. 2016;44(6):1086–105.
8. Staudenmayer J, Zhu W, Catellier DJ. Statistical considerations in the analysis of accelerometry-based activity monitor data. Medicine & Science in Sports & Exercise. 2012;44(1 Suppl 1):S61–7.
9. Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet. 1986;1(8476):301–10.
10. Giavarina D. Understanding bland altman analysis. Biochemia Medica. 2015;25(2):141–51.
11. Neugebauer JM, Collins KH, Hawkins DA. Ground reaction force estimates from ActiGraph GT3X+ hip accelerations. PLoS One. 2014/06/11. 2014;9(6):e99023.